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Abstract. The 2008 NRDA conference introduced the Numerical INJection 
Analysis project (NINJA), a new collaborative effort between the numerical 
relativity community and the data analysis community. NINJA focuses on 
modeling and searching for gravitational wave signatures from the coalescence 
of binary system of compaet objects. We review the scope of this collaboration 
and the components of the first NINJA project, where numerical relativity groups 
shared waveforms and data analysis teams applied various techniques to detect 
them when embedded in colored Gaussian noise. 



1. Introduction 

The coalescence of binary systems of c;onipac;t objects is one of the most promising 
sources of gravitational wave radiation for ground-based detectors. The past few years 
have witnessed the successful construction and operation of a world-wide network 
of large interferometric gravitational wave detectors, with the three LIGO detectors 
in the United States [1], Virgo in Italy [2], TAMA in Japan [3] and GEO600 in 
Germany [4] . A number of searches for unmodelled bursts and binary coalescences are 
in progress and several results have already been reported [5-22]. In parallel, recent 
breakthroughs in numerical relativity made it possible to successfully simulate the 
merger phase of Binary Black Hole (BBH) coalescences. Following the pioneer work 
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of Pretorius [23] , and the independent results from the Goddard [24] and Brownsville 
(now at RIT) [25] groups, several numerical relativity groups around the world have 
successfully evolved various BBH configurations. This has lead to important new 
physical insights, such as the prediction of large recoil velocities due to asymmetric 
emission of gravitational radiation during the merger [26-42] and the prediction of the 
parameters of the remnant for a wide class of initial states [43-54]. A review of the 
current status of BBH simulations is available in this volume [55] . 

Numerical relativity now provides a complete model of the coalescence waveform, 
inclusive of the merger. However, the numerical relativity results have not yet been 
synthesized into an analytical model over a broad range of mass ratios, spin, and 
eccentricity. Furthermore, despite significant progress, there is not yet a complete 
description of how post-Newtonian and numerical simulations are to be matched with 
each other, over the full parameter space. Most searches for gravitational waves from 
BBH mergers have so far relied on the post-Newtonian analytical description of the 
expected gravitational wave signal, valid only when the black holes are sufficiently far 
apart. A verification of their robustness is urgent, and, more generally, it is important 
to quantify the performance of the data analysis pipelines for detection and parameter 
estimation, to set astrophysical upper limits, and to aid in follow-up studies and the 
interpretation of potential detection candidates. 

The Numerical INJection Analysis project (NINJA), open to all scientists 
interested in numerical simulations and gravitational wave data analysis, was started in 
the spring of 2008 to address these issues and form a new, close collaboration between 
the two communities. The first NINJA project benefited from the contribution of 10 
numerical relativity groups and 9 data analysis teams, for a total of 76 participants 
from 30 institutions. In a step towards the incorporation of numerical relativity 
waveforms in gravitational wave data analysis, various data analysis algorithms 
analyzed BBH coalescence waveforms buried in simulated Gaussian noise at the design 
sensitivity of initial LIGO and Virgo. The first NINJA project only includes BBH 
simulations, but it is expected that future NINJA analyses will be expanded to include, 
for example, binary neutron star and supernovae simulations. 

The first NINJA project was a learning experience which brought the two 
communities closer. To facilitate this, minimal constraints were imposed: each 
numerical relativity group chose which waveforms to share and each data analysis 
team chose which methods and results to contribute. We learned about technical and 
conceptual issues, and how to address them. However, due to the small statistics of 
simulations and a lack of systematic studies, comparisons and conclusions drawn from 
the first NINJA project are limited and should be handled with care, as in most cases 
they are only the first steps towards fully understanding the sensitivity of data-analysis 
pipelines to black hole signals. 

This proceedings paper provides a highlight overview of the NINJA collaboration 
and its first project. For details on the waveforms, on the data analysis techniques 
and the results, we refer to [56]. 

2. Numerical Relativity Waveforms 

Ten numerical relativity teams contributed BBH coalescence waveforms that are 

solutions to Einstein's equations, with no restrictions on their morphology or accuracy. 
Each team submitted a maximum of two waveforms, or up to five waveforms if they 
were part of a one-parameter family, following the format specifications in [57]. 
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The NINJA waveforms cover a variety of physical and numerical parameters. 
Most simulations model low-ccccntricity inspiral, with mass ratio q = M1/M2 from 
1 to 4, and several spin configurations; the initial frequency of the £ = m = 2 
mode ranges from 0.033/M to 0.203/M, where M is the sum of the initial black-hole 
masses, and the waveform length varies between a few lOOM and over 4000M. The 
different contributions are: BAMJiHB [58-62] and BAM_FAU [53,58,59,62], using 
the BAM code, the AEI/LSU code CCATIE [29, 37, 46, 63, 64], the Goddard Space Flight 
Center's code Hahndol [65, 66], the RIT code LazEv [25, 40, 67], Ulrich Spcrhakc's code 
Lean [50, 68, 69], the Georgia Tech/Penn State code MayaKranc [51, 70], the Princeton 
University code [23,71-73], the Cornell/Caltech collaboration code SpEC [74-77], and 
the University of Illinois at Urbana-Champaign code [78]. 

The numerical codes follow one of two approaches to solving the Einstein 
equations: the generalized harmonic formulation [23], or the moving-puncture 
approach [24,25]. All of the results make the simplifying assumption of conformal 
flatness for the spatial metric of the initial slice, which leads to some spurious 
gravitational radiation in the initial data, and attempt to model non-eccentric 
inspirals, with the exception of PU T52W and MayaKranc c02 that target eccentric 
inspirals. To estimate the gravitational-wave signal at a finite distance from the 
source, the SpEC and CCATIE contributions use the Zerilli-Moncrief perturbative 
formalism [79 81], while all others use the Newman- Penrose curvature scalar t/j^ [82]. 
Details on the implementations within particular codes can be found, for instance in 
refs. [37, 58, 68, 83] ; a detailed review of similarities and differences between NINJA 
waveforms, codes and numerical methods, is available in rcf. [56]. 

Waveforms have been contributed in the form of spherical harmonic modes ~'^Yi^ 
of spin- weight —2 of the radiation field at large distance from the source, as specified 
in [57]. In the Transverse- Traceless (TT) gauge, the spatial components htj are: 

hij=Aijy+0{r-^), (1) 

where M is the total mass of the system, r is the distance from the source, and Aij is 
a time dependent TT tensor. In the TT gauge, hij has two independent polarizations 
/i+ and hx, and the complex function — ihx is decomposed as: 

M °° ^ 

h+-ihx = —Y^ 51 ^^"^ (*) ~^"^irn{t, (t>) ■ (2) 

The expansion parameters Him are complex functions of the retarded time t — r and, 
if r is the radius of extraction of the wave. Him is function of t only. The angles l 
and (j) are respectively the polar and azimuthal angles in a suitable coordinate system 
centered on the source. Computation of the strain from the Zerilli-Moncrief odd- 
and even-parity (Q;^, Qtm) multipoles of the metric perturbation requires one time 
integration [37,80]; computation of the strain from the Newman- Penrose curvature 
scalar ifj^ requires two time integrations, with a proper choice of the constants of 
integration, and attention to artifacts resulting from the finite extraction radii. 

While no attempt was made, in this first project, to restrict waveform parameters, 
future exchanges will be coordinated to address questions such as how many modes 
are needed for complicated waveforms: should a fixed i cutoff be imposed or should 
all modes with a minimum percentage fraction of the total energy be used? What 
is the impact of this on the total systematics as a function of mass and the detector 
noise curve? We will also explore hybridization with post-Newtonian waveforms, to 
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avoid abrupt startup of the waveforms, and exploit overlaps and complementarities 
between teams and codes, for a better coverage of the parameter space. 

3. The NINJA Data Set 

Following the format outlined in [57] , numerical relativity groups provided one ASCII 
data file for each mode {£, m) with appreciable contribution to the waveform. Each 
data file consists of three columns: time in units of the total mass, and the real and 
imaginary parts of the mode coefficients 

The NINJA data set simulated signals as seen by the gravitational wave detectors, 
which also depend on the total mass of the binary system, its distance, and its 
orientation with respect to a particular detector. The data production starts from 
a Monte Carlo population of binary black hole systems, with the source distance 
logarithmically distributed in 50 — 500 Mpc. Orientations are uniformly distributed 
and masses are uniform in 20 — 35OM0. Since the total mass M scales both time 
and amplitude, each waveform can be scaled to an arbitrary value of total mass, with 
the requirement that the initial frequency of the dominant £ = m = 2 mode is below 
30 Hz, to avoid startup transient artifacts. This is an appropriate threshold for the 
sensitivity curve of the LIGO and Virgo initial detectors, and sets a lower limit on 
the value of M for each waveform: the longer the waveform, the lower the value of 
M it can be injected at. No restrictions were placed on other simulation parameters, 
such as spin, mass-ratio and eccentricity, as these were determined by the numerical 
waveform parameters. 

Given Hi^nit), the total mass, the distance to the source, the angles (t, ^), /i+,x 
(from eq. 2), and the detector response functions -F+,x the observed strain is: 

h{t) = h+{t)F+{Q, V) + /ix {t)F^ (e, V) , (3) 

where (6, are sky angles in the detector frame and ij) is the polarization angle. The 
time t is in units of seconds. 

Once signals for the chosen populations are produced as time series at 4096 Hz 
sampling rate, typical for LIGO and Virgo searches, they are added to stationary 
Gaussian noise colored with the initial LIGO and Virgo design power spectral densities. 
A data set was produced for each of the three LIGO interferometers and for the Virgo 
interferometer, with approximately 12 signals from each contributing group. The 
resulting NINJA data set consisted of 126 signals injected in a little over 30 hours; the 
time interval betwen adjacent signals is a random number between 700 ± 100 sec. The 
distribution of signal-to- noise ratios is detector-dependent; Fig. 1 shows SNR versus 
total mass at the LIGO Hanford 4km detector. 

4. Data Analysis 

Nine data analysis teams contributed a variety of techniques to the NINJA project. 
Participants were provided with the data set described in section 3 and a list of 
injected signals parameters, but did not have access to the raw waveforms. The two 
fundamental data analysis goals are detection and parameter estimation. Detection 
has been pursued in the NINJA project with both matched-filtering and un-modeled, 
excess-power techniques. Bayesian model and parameter estimation techniques were 
also applied to the NINJA data set. In this section we provide an overview of the 
techniques adopted; for more details and results, see [56]. 
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Figure 1. Distribution of signal-to-noisc ratio versus total mass for the NINJA 
data set. The SNR is detector-dependent; in this plot, it is computed at the 
LIGO-Hanford 4km detector. 



Matched filter, modeled searches 

The LIGO Scientific Collaboration (LSC) and the Virgo Collaboration have 
implemented a matched filter analysis pipeline in previous searches for gravitational 
wave signatures from compact binary coalescences [8,84]. For each waveform model, 
a bank of templates covers the parameter space so that the fractional loss in signal- 
to-noise ratio (SNR) between any signal and the nearest template is less than 3%. 
Data from each detector is match filtered against this bank [85, 86] and a trigger is 
produced whenever the SNR exceeds a threshold of 5.5. Triggers which do not have 
coincident parameters in two or more detectors are discarded [87, 88] ; the remaining 
triggers are ranked by the square-sum of the SNRs of the triggers in the coincidence. 
For inspiral matched filter analyses, coincident inspiral triggers are subject to a second 
filtering stage, in which signal-based vetoes are also applied to separate true signals 
from noise fluctuations [89, 90]. A threshold is applied to an effective SNR pcff , which 
combines the matched filter SNR and the value of the signal-based veto [89]. 

Six groups contributed matched filter searches to the NINJA project; the results 
can be divided into three categories based on the waveform templates used. Within 
these categories, different parameter choices were made to explore how the pipeline 
can detect numerical relativity simulations. 

4.. 1.1. Inspiral only waveforms from the Post- Newtonian expansion. The standard 
template for the LSC-Virgo search [5-8, 10] is based on the non-spinning post- 
Newtonian inspiral waveforms, calculated directly in the Fourier domain through the 
stationary phase approximation [85,91]. These waveforms, referred to as SPA or 
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TaylorF2, are parameterized by the binary's component masses mi and m2 [92-94], 
or cquivalcntly by total mass M = mi +m2 and symmetric mass ratio rj = mim2/M^. 
The amphtude evolution is modelled to leading order and the phase evolution is 
modelled to a specified post-Newtonian order. Formally, the TaylorF2 waveform 
can be extended to arbitrary frequencies. However, at higher frequencies, the post- 
Newtonian expansion does not accurately model the physics and the waveforms are 
terminated at a cutoff frequency fc- In the LSC- Virgo analyses, this is chosen to be the 
innermost stable circular orbit (ISCO) frequency for a test mass in a Schwarzschild 
spacetime. Details of a search of the NINJA data using TaylorF2 waveforms are 
available in this volume [95] . 



4-1-2. Ring down- only waveforms from black hole perturbation theory. Ringdown 
searches use a two-parameter template bank parameterized by frequency and quality 
factor, constructed to cover the desired range of mass and spin [86]. The LSC ringdown 
matched filter algorithm [88] was used for the NINJA analysis, using a bank with 
frequency between 50 Hz and 2 kHz and quality factor between 2 and 20. 

Assuming that the dominant ringdown mode is the fundamental mode, £ = m = 2, 
the frequency and quality factor of the waveform can be expressed in terms of the 
mass M and dimensionless spin factor a of the black hole using Echeverria's fit [96] 
to Leaver's numerical calculations [97]: 

f=l^b-^-^Hl-a)^ (4) 

= 2(l-a)-w. (5) 

These equations can be inverted to calculate the mass M and spin a from the template 
parameters of a given trigger. 



4.. 1.3. Phenomenological/ analytical models for the full coalescence. Post-Newtonian 
theory is valid when the black holes are sufficiently separated, but becomes unreliable 
as the velocity of the black holes increase in their final orbits before merger and the 
non-perturbative information contained in numerical simulations becomes important. 
A successful approach has been to combine analytical (PN) and numerical (NR) results 
into full waveform templates. Three different families of waveforms have been used 
in analyzing the NINJA data: Extended TaylorF2, EOBNR and phenomenological 
waveforms. 

Recent comparisons to numerical relativity waveforms have shown that extending 
the TaylorF2 waveforms to higher frequencies improves their sensitivity at higher 
masses [98,99]. An analysis using these extended waveforms was performed on the 
NINJA data. Specifically, the cutoff frequency /c was increased from ISCO to the 
effective ringdown (ERD) frequency, obtained by comparing post-Newtonian models to 
numerical waveforms [98]. This increases the sensitivity of the search as the extended 
templates can detect some of the power contained in the late inspiral or early merger 
part of the signal. 

Two studies of the match between numerical waveforms [98,99] and TaylorF2 
templates suggest that the search efficiency can be improved by extending the range 
of r] outside the physical range of r; < 0.25, and terminating the waveforms at a 
weighted ringdown (WRD) frequency, between ISCO and ERD [99]. By extending 
the template bank to cover all points with r/ < 1, its size was approximately doubled. 
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An analysis of the NINJA data with the extended template bank showed an increse 
in the recovered SNR for some signals. 

Combining results from PN and perturbation theory, the EOB model [100, 101] 
predicted inspiral, merger and ringdown waveforms. The non-spinning EOB model has 
been further improved by calibrating it to NR results, achieving rather high matching 
perfomances without maximizing on binary parameters, but only on initial phase and 
time of arrival [72, 98, 102-106] . These waveforms, called EOBNR, were also used in 
analyzing the NINJA data. 

The phenomenological approach [107, 108] matches PN and NR waveforms in a 
regime where both are sufficiently accurate, and fits the resulting hybrid waveform to a 
phenomenological model in frequency domain. This procedure has been carried out for 
non-spinning black holes; each waveform is parametrized by the physical parameters 
of the system, i.e. the masses mi and m2 of the black holes. Results of a search of 
the NINJA data with phenomenological templates are given in [109]. 

4-2. Unmodeled searches 

Burst searches are designed to detect gravitational wave transient signatures with 
minimal assumptions on their origin and waveform. They do not use templates 
and instead target excesses of power in the time-frequency plane. Since they do not 
assume a template, and they target short transients, burst searches are suited for the 
detection of the merger phase of the coalescence, and have the potential to probe a 
large parameter space, inclusive of spin and ellipticity, at no additional computational 
cost. Two such algorithms, developed by the LSC and Virgo, analyzed the NINJA 
data: Q-pipeline and HHT. 

The Q-pipeline [110-112] is a multi-resolution time- frequency search for 
statistically significant excess signal energy, equivalent to a templated matched filter 
search for sinusoidal Gaussians in whitened data. The template bank covers a 
finite region in time, frequency, and quality factor such that mismatch between any 
sinusoidal Gaussian in this signal space and the nearest basis function does not 
exceed a maxiimiin of 20% in energy. For the purpose of the NINJA project, the 
Q-pipeline analysis focused on the detection efficiency at the single detector, with an 
SNR threshold of 5.5, as in the matched filter searches. 

The Hilbert- Huang Transform (HHT) [113 115] is an adaptive algorithm that 
decomposes the data into Intrinsic Mode Functions (IMF), each representing a unique 
locally monochromatic frequency scale of the data. The original data is recovered 
with a sum over all IMFs. The Hilbert transform, applied to each IMF, provides 
instantaneous frequency and amplitude as a function of time, thus providing high 
time-frequency resolution for signal detection. Details on how this pipeline was tested 
on NINJA data are available in this volume, in ref. [116]. 

4-3. Comparison of searches 

The analysis techniques applied to NINJA data have comparable performance, if the 
same SNR threshold is used. However, the simulation statistics in this first NINJA 
project were quite limited, and the SNR of injected signals was relatively large, as 
shown in Fig. 1. Gonsequently, no strong conclusions can be drawn from these results 
except that unmodeled algorithms and matched filtering to ringdown or to the "wrong" 
templates can still identify loud signals in Gaussian noise, although not necessarily 
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with the correct parameters, as they preferentially detect whichever portion of the 
signal is in the detectors' most sensitive band. For the lowest injected masses, this 
corresponds to the inspiral phase and the detected frequency is close to fiscOi but as 
the injected mass increases, the ringdown dominates, as expected. For a quantitative 
assessment of these comparison, with details, we refer to [56]. 



4.^. Bayesian parameter estimation and model selection 

Bayesian inference is a powerful means of extracting information from observational 
data; although computationally too expensive for a search, it can be very useful for 
parameter estimation, in the study of candidates identified by search pipelines that 
are based on a pre-determined template bank. 

Two (iifferent approaches were taken to make inferences on the NINJA data: (i) 
the estimation of the parameters of the signal assumed to be present in the data, and 
(ii) the calculation of confidence in the presence of the signal, quantified by the odds 
ratio between the signal and noise models of the data. 

Both approaches require the calculation of the posterior probability-density 
function (PDF) on the parameter space of the signal, given the data d: 

. . / ~ ~ . 2 

Pi9)pid\9) 



d{f)~h{f;0) 



Snif) 



■df\ (6) 



in the presence of Gaussian noise with power spectral density Sn{f), where p{6) is the 
prior probability density of the parameters 9 and h{9) is the model used to describe 
the signal. 

A Markov-Chain Monte-Carlo (MCMC) method [117, 118] was used to coherently 
analyze data from multiple detectors and evaluate the posterior PDFs. This technique 
stochastically samples a 12-dimensional parameter space for the best match to the 
data. The accuracy of the parameter estimation is assessed with a comparison of 
likelihoods at difi^erent points in parameter space. Inspiral waveforms at 1.5 post- 
Newtonian order were used in the NINJA analysis, with good parameter estimation 
for injections with relatively low total mass, where the inspiral contained a significant 
fraction of the total signal-to-noise ratio. For high-mass injections, matching the 
merger and ringdown portions of the waveform to inspiral templates resulted in poor 
parameter estimation. More details on the MCMC code are available elsewhere in this 
volume [119]. 

Bayesian model selection, based on the nested samphng [120] Monte Carlo 
technique, was used to measure the confidence of a detection with a given template 
family. In this approach, the likelihood of signal and noise model are marginalized 
over a range of parameters, with their prior. 

The ratio of probabilities, or posterior odds ratio, of the two models is the product 
of the prior odds ratio and the Bayes factor, computed from the data (eq. 6) and 
corresponding to the level of confidence in the detection [121,122]. Bayes factors 
were computed for the injections in the NINJA data set with a model that inc;luded 
the coherent response of all four available interferometers and two non-spinning 
approximants: the standard (2 PN) TaylorF2 family, with inspiral truncated at 6M, in 
50— 150Mq, and the phenomenological inspiral-merger-ringdown templates described 
in [107] in 50 — 475 M©, with uniform priors for all parameters. For all the injections, 
the phenomenological approximants return a Bayes factor significantly larger than 
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TaylorF2 approximants, which do not contain the merger and ring-down portion of 

the coalescence. 

The algorithm also yields the maximum likelihood estimate of the recovered 
parameters, with a measure of the statistical errors. The resulting total mass is in 
most cases underestimated, but the sky location is on average fairly well determined, 
as it depends on the time of arrival of the signals at different instrument sites. This is 
currently under more careful investigation, for more detail please see, in this volume, 
ref. [123]. 

5. Conclusions 

The NINJA project was conceived a first step towards a long term collaboration 
between numerical relativists and data analysts in an attempt to use numerical 
waveforms to enhance searches for gravitational waves and bridge the gap between 
gravitational wave observations and astrophysics. The work is unique in that it focused 
on running existing gravitational wave search algorithms on waveforms from numerical 
simulations. To facdlitatc this first collaboration and encourage broad participation, 
the scope was limited, with a relatively small number of waveforms exchanged and 
limited coordination between groups. As a result, the statistics of simulations was 
too low to draw robust c;onc;lusions [56]. NINJA has proven to be extremely valuable 
at framing the questions which need to be answered for an effective use of numerical 
relativity waveforms in data analysis: how many modes are needed for complicated 
waveforms? What accuracy is appropriate for each analysis? How important is 
hybridization with post-Newtonian waveforms? How can we best cover the parameter 
space? Which approximate waveform families are suitable for the detection of nature 's 
gravitational-wave signals, as modeled by numerical-relativity waveforms, without 
significant loss of signal-to-noise ratio? Which approximate families, if any, are 
adequate for accurate parameter estimation? Future NINJA analyses will build upon 
this work to address these questions. 
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